Phase diagram of Hertzian spheres 
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We report the phase diagram of interpenetrating Hertzian spheres. The Hertz potential is purely 
repulsive, bounded at zero separation, and decreases monotonically as a power law with exponent 
5/2, vanishing at the overlapping threshold. This simple functional describes the elastic interaction 
of weakly deformable bodies and, therefore, it is a reliable physical model of soft macromolecules, 
like star polymers and globular micelles. Using thermodynamic integration and extensive Monte 
Carlo simulations, we computed accurate free energies of the fluid phase and a large number of 
crystal structures. For this, we defined a general primitive unit cell that allows for the simulation 
of any lattice. We find multiple re-entrant melting and first-order transitions between crystals with 
cubic, trigonal, tetragonal and hexagonal symmetries. 



I. INTRODUCTION 

Many soft materials consist of building blocks that 
are themselves "soft", meaning that, under experimen- 
tal conditions, they can interpenetrate or deform appre- 
ciably due to the forces acting between them. Linear 
polymers provide an example!^ at moderate free-energy 
cost, the centers of mass of two distinct polymers can be 
made to superimpose. These macromolecules are there- 
fore fully interpenetrable. Other particles do not inter- 
penetrate but are deformable, like polymer micelles: they 
deform elastically at high densities, thus exploring con- 
figurations that are forbidden to hard-particle systems. 

Systems of soft particles are interesting because crys- 
tallization in simple and colloidal systems is usually at- 
tributed to excluded-volume effects. Systems with soft- 
repulsive, short-range interactions are therefore expected 
to have a qualitatively different freezing behavior than 
those with hard particles. Indeed, one of the striking 
features of systems of particles with a bounded repulsive 
interaction is that they have a maximum melting temper- 
ature. Above it, the solid phase is unstable. This phe- 
nomenon can be understood as follows: consider a sys- 
tem with an interparticle potential that has its maximum 
energy e at zero separation and decreases monotonically 
to zero within a distance a. When the thermal energy 
ksT ^ e, interparticle interactions become negligible 
and the system approaches ideal gas behavior. However, 
when ksT <^ e particles can only overlap at a significant 
energetic cost. In the limit ks T/e 0, the system will 
undergo a hard-sphere freezing transition. In fact, some 
30 years ago, Stillinger already showed that the so-called 
Gaussian Core Model (GCM) has a maximum melting 
point M oreover, Stillinger, and subsequently other au- 
thorSjl^ showed that the GCM also exhibits re-entrant 
melting upon compression, and a transition from face- to 
body-centered cubic crystal structures. Such a rich phase 
behavior is not limited to the repulsive Gaussian poten- 
tial, but has been shown to be common to many effective. 



soft potentials that were designed to reproduce the phase 
behavior and dynamics of microgelSj'^star-polyelectrolyte 
solution^ and star-polymer solutions .l^l 

However, not all soft-repulsive potentials give rise 
to re-entrant melting. The penetrable sphere model 
(PSM),'i2l a square-shoulder potential, is the simplest 
model that leads to freezing at any arbitrary temperature 
and to the appearance of cluster crystals (multiple occu- 
pation of lattice sites). Multiple occupancy in the PSM 
arises because, as soon as two particles interpenetrate, 
no further energetic cost exists for full overlapping. As a 
consequence, particles can sit on top of each other thus 
reducing the total number of overlaps and the interaction 
energy. Density functional theory and molecular simula- 
tions have confirmed the occurrence of cluster crystals for 
the generalized exponential model.lSl It is expected that 
certain dendrimers could form such crystalSjI^^l but thus 
far experiments are lacking (see, however, ref. I13p . 

It is apparent from the discussion above that the 
shape of the pair potential determines when the re- 
entrant melting and clustering scenarios manifest in the 
phase diagram. Indeed, Likos et al^ established a 
simple criterion to distinguish between the two scenar- 
ios for bounded repulsive potentials: re-entrant melting 
happens for bounded potentials with a positive definite 
Fourier transform. Otherwise, clustering and freezing are 
expected to occur at all temperatures. 

Although the type of scenario, either re-entrant melt- 
ing or clustering, can be predicted, there is no method 
that guarantees an a priori prediction of all stable crys- 
tal structures given the bounded repulsive potential. The 
usual method involves selecting a small number of can- 
didate structures, for which the ground-state energy or 
finite-temperature free energies are computed. Among 
the candidates, the stable structure at any particular 
density and temperature is the one with the lowest (free) 
energy. This works well for typical pair potentials, for 
which one assumes that unusual crystals are unlikely to 
be stable structures. However, soft potentials are known 
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to induce multiple crystals)^ including the less typical 
non-cub ic lattic es, and exotic structures like the A15 and 
diamoncPSmEII (for an extreme example, see ref. [TS). 

There exist powerful approaches to search for stable 
crystal structures, such as those based on metadynam- 
icpSl and genetic algorithmspH However, metadynamics 
can be very expensive, in particular when trying to es- 
cape deep, local free-energy minima. In a genetic algo- 
rithm, an initially random population of lattices is mod- 
ified using evolutionary rules according to a fitness func- 
tion: high fitness corresponds to low (free) energy. The 
structure with the highest fitness in the population after 
breeding a few generations is assumed to be the most sta- 
ble lattice. Although genetic algorithms greatly alleviate 
the problem of searching for candidate crystal structures, 
they still suffer from the same numerical bottleneck of the 
usual method: the expensive computation of free ener- 
gies. For this reason, numerical searches for stable crystal 
structures are often based on energies, rather than free 
energies. To improve upon this approach, one can ap- 
proximate the free energy by that of a harmonic crystal 
or by using density- functional theory. However, to im- 
prove accuracy, one has to turn to computer simulation 
methods that generate free energies that are "exact" to 
within statistical accuracy. Accurate free-energy calcu- 
lations can be crucial in soft systems where, specially 
at high densities, crystals with very similar energies can 
compete for stability. 

In this paper we focus on the phase behavior of 
Hertzian spheres. The Hertz potential describes the 
change in elastic energy of two deformable objects when 
subjected to an axial compression.^^ The potential has 
the following form: 

Vir)^[^^'-f^''' :<^, (1) 

where a and e set the length and energy scales and r is the 
distance between the centers of the undeformed spheres. 
In the limit ksT/e — > the hard-sphere model is recov- 
ered. The Hertz model was developed to describe the in- 
teraction between elastic spheres that are deformed only 
slightly. Clearly, for very large compressions or at high 
densities (well above the overlap concentration) the as- 
sumptions underlying the original Hertz model no longer 
apply. Here we use the Hertz potential as a simple repre- 
sentation of soft particles that is bounded (i.e. it remains 
finite at r = 0), has a finite range (viz. a) and a positive, 
yet sharply decaying Fourier transform. Also, the Hertz 
potential is computationally cheap and thus allows us to 
perform extensive simulations for the calculation of free 
energies and coexistence points. We have not attempted 
to map the present model onto a specific soft-colloid or 
star-polymer system.EEl just note that a judicious de- 
sign of star polymers or dendrimerdf^I makes it possible 
to "engineer" the interaction between these particles to 
a considerable degree. 

Using standard thermodynamic integration, we calcu- 
lated free energies of the fluid phase and of candidate 




FIG. 1: The general primitive cell and vectors pi, p2 and ps. 

crystal structures. As candidates, we considered all Bra- 
vais lattices and the hexagonal close-packed, diamond, 
and A15 structures. Our results show that the phase 
diagram of the Hertz model exhibits multiple re-entrant 
melting and polymorphic transitions between a number 
of crystals. In addition, we find that the re-entrant fluid 
shows unusual diffusivity curves, and speculate about the 
nature of this fluid at high density and low temperature. 



II. GENERAL PRIMITIVE UNIT CELL 
A. General framework for Bravais lattices 

All 14 Bravais lattices can be generated from a general 
primitive unit celP^^ defined by three independent vec- 
tors pi, P2 and P3 G M^. Once these primitive vectors 
are defined for a particular crystal, the location of every 
particle in the crystal follows from 

r = nipi -I- n2P2 + riapa -I- ro, (2) 

where ni,n2 and 713 G Z, and Tq is an arbitrary vector. 
For simplicity, we can fix one of the particles at the origin 
of the coordinate system, and thus Tq = (0, 0, 0). 

There are many ways of defining primitive vectors for 
each Bravais lattice. Here we conveniently choose the 
general set of vectors depicted in Fig. [T] where pi is par- 
allel to the X axis, and pi and P2 lie on the xy plane of 
the cartesian coordinate system. Algebraically, 

-Pi = (1,0,0), (3) 
a 

-P2 = ^ (cos 7, sin 7,0), (4) 
a a 

1 c 

-P3 = - (cos COS a' sin /3, sin a' sin/?), (5) 
a a 

where a, b, and c are the lengths of pi, P2 and P3, re- 
spectively, and a, (3 and 7 denote the angles formed by 
each pair of primitive vectors, as indicated in Fig.[l] The 
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TABLE L Parameters characterizing the 14 Bravais lattices according to Eqs. |2)8[ Round and square brackets indicate free 
and dependent parameters, respectively. 



Bravais lattice 


b/a 


c/a 


Qf/vr 


f3/n 




V 


oimpie cuDic i^ov.^ ) 


1 


1 


1 /9 


1 /9 


1 /9 
l/Z 


1 

i. 


Face-centered cubic (FCC) 


,/9 /9 


,/9 /9 


1 


i/4 


1 /A 
1/4 


1 /A 
1/4 


Body-centered cubic (BCC) 


1 
1 


\J 6/ Z 


U.OU4 


U.OU4 


1 /9 
1/^ 


1 /9 
i/^ 


Hexagonal (H) 


-I 
1 




1 /9 


1 /9 
1/2 


1/3 


[V 3/ 4J 




\ 


1 /2 


1 /2 

i/Z 


1 /2 


1 /2 

i/Z 


1 /2 

1/ z 


Body-centered tetragonal (BCT) 


1 


(3/4) 


\0 2681 


10 2681 


1/2 


11/41 


Rhombohedral or trigonal (R) 


1 


1 


(1/3) 


[1/3] 


[1/3] 


[^2/2] 


Simple orthorhombic (SO) 


(V2/2) 


(^/3/2) 


1/2 


1/2 


1/2 


[0.612] 


Base-centered orthorhombic (BaCO) 


(1) 


iV2/2) 


1/2 


1/2 


[1/4] 


[1/2] 


Body-centered orthorhombic (BCO) 


(^/2) 


(4/5) 


[0.285] 


[0.285] 


1/2 


[0.265] 


Face-centered orthorhombic (FCO) 


(3/5) 


(^/3/2) 


[0.340] 


[0.304] 


[0.186] 


[0.235] 


Simple monoclinic (SM) 


(V2/2) 


(V3/2) 


(1/3) 


1/2 


1/2 


[0.530] 


Base-centered monoclinic (BaCM) 


1 


(^/2/2) 


(1/3) 


1/3 


(1/4) 


[0.420] 


Triclinic (T) 


(V2/2) 


(V3/2) 


(1/3) 


(0.4) 


(2/3) 


[0.306] 



TABLE II: Basis vectors for the non-Bravais lattices considered in this work. 



Non-Bravais lattice Reference Bravais lattice nt hi /a 



Diamond (D) FCC 1 

Hexagonal close-packed (HCP) Simple hexagonal (SH) (H with c/a = 1) 1 
A15 SC 7 (ii|)(|,0,i)(iO,|)(il,0) 



^ 1 1 1 ^ 

^4 ' 2 ' 2 J 
k 2 ' 3 ' 2 J 



;i,i,o) (o,i,i) (0,1,1) 



angle a' is such as 



cos a 



sma 



cos a — cos f3 cos 7 
sin /? sin 7 ' 



\/— cos^ a — cos^ /3 -I- sin^ 7 + 2 cos a cos /? cos 7 
sin /3 sin 7 



(6) 



(7) 



The volume V of the primitive cell is |pi • (p2 x pa) 
and, therefore. 



V be . , 

= sma sinpsm7. 



(8) 



From eq. [8] we can see that the length a sets the spe- 
cific volume of the crystal. Without loss of generality, 
we choose a to be the longest among the three primitive 
vectors. Hence, according to eqns.[2][7] the 14 Bravais lat- 
tices can be described with five parameters: < 6/a < 1, 
< c/a < 1, < a < tt/2, < /? < 7r/2 and < 7 < tt. 
Table |T] shows the specific values of the parameters for 
the 14 Bravais lattices. Depending on the particular lat- 
tice, any of the five parameters can be free, fixed or de- 
pendent. A free parameter, shown in round brackets in 
Table |T] can be changed at will to generate different spe- 
cific structures belonging to the same lattice, provided 
that — cos^ a — cos^ f3 + sin^ 7 -I- 2 cos a cos /? cos 7 > 0. A 
dependent parameter, shown in square brackets, has val- 
ues that depend on those assigned to the free parameters. 



according to the expressions in Table III By definition, 
some lattices are a subset of the tetragonal, trigonal, or- 
thorhombic, monoclinic or triclinic, as indicated in Ta- 
ble lIVl 



B. Extension to non-Bravais lattices 



Bravais lattices have 1 particle in their primitive unit 
cell. The primitive unit cell of non-Bravais crystals con- 
tains rib additional particles that are characterized by a 
set of translation or basis vectors b^, where i runs from 
to rif,. Because of translational symmetry, one of the 
basis vectors can be chosen arbitrarily, and this is why 
in eq. [2] we defined bp = rg = (0,0,0). A non-Bravais 
crystal is generated by translation of its reference Bravais 
lattice using the set of basis vectors 



(9) 



4=0 



Table llTl shows the reference lattice and basis vectors cor- 
responding to the three non-Bravais crystals considered 
in this work. 
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TABLE III: Lattice-dependent parameters. 



TABLE IV: Parameters leading to equivalent lattices. 



Bravais lattice 



Dependent parameters 



BCT, BCO 
R 

ECO 
ECO 

BaCO, ECO 



Q = arccos 



a = (3 = arccos [^2^ j 
f3 — 'y — a 



/3 = arccos |^2^j 
7 = arccos 



C. General periodic boundary conditions 



According to the framework described in sections |II A| 
and |IIB[ one can generate any crystal structure as a pe- 
riodic repetition of the primitive unit cell along its three 
primitive vectors. Equivalently, one can also construct a 
unit cell that consists of an arbitrary number of primi- 
tive cells, and replicate the unit cell along its lattice vec- 
tors to obtain the same periodic structure. Therefore, we 
can define a non-primitive unit cell that contains a suffi- 
ciently large number of particles, and use general periodic 
boundary conditions to simulate the bulk crystal. The 
lattice vectors that define this unit cell are Ui — nipi, 
U2 = n2P2 and U3 = 713P3, where ni, ri2 and are 
the number of primitive cells contained in the unit cell in 
each of the directions of the lattice vectors. 

The nearest image of a particle in a unit cell can be 
found using the reciprocal lattice vectors Vi, V2 and V3, 
which are defined such as • Vj = 6ij . Using the general 



primitive vectors of eqns. |3][5j the corresponding recipro- 
cal vectors can be written as 



a vi = 



av2 



av3 = 



U2 X U3 



(U2 X U3) • Ui 

— cos 7 cos a' sin P cos 7 — cos (3 sin 7 
sin a' sin (3 sin 7 
1 /„ 1 —cos a' 



1, 



sin7 

U3 X Ui 



(U3 X Ui) • U2 b/a \ ' sin 7' sin a' sin 7 



0, 



Ui X U2 



(ui X U2) • U3 c/a 



0,0,^ 



1 



sin a' sin (3 



, (10) 

,(11) 
(12) 



Following eq. [2] the location r„ of a particle inside the 
unit cell with an image at position r can be calculated as 



where 



r - rTiiUi - m2U2 - m3U3, 



= {r • Vj} for i = 1,2,3 



(13) 



(14) 



and { • • • } denotes the nearest integer. Note that, by 
definition, {r,; • v^} = 0. 



Bravais lattice 


Parameters and equivalent lattice 


ST 


c/a = 1 SC 


BCT 


c/a = x/3/2 ^ BCC 


R 


SC 


SO 


b/a = 1 or c/a = 1 or 6/a = c/a i-+ ST 




b/a = c/a = 1 i-> SC 


BaCO 


b/a = c/a = V2/2 SC 


BCO 


6/a = 1 BCT 




b/a = 1 and c/a = V3/2 ^ BCC 


ECO 


b/a = c/a = V2/2 ^ FCC 


SM 


a^n/2^ SO 


BaCM 


c/a = 1 and a — "f — 1/3 t-^ R with a — n/3 


T 


parameters of Table Pll t-^ any of the lattices 



III. THERMODYNAMIC INTEGRATION 

In order to construct the phase diagram, we first per- 
formed a rough test of the stability of the fluid phase 
by running short NVT Monte Carlo simulations at 100 
points equally distributed in the T — p space, T/e 
within the interval [0 — 0.01] and pa^ within [0 — 10]. 
This allowed us to map approximately the solid region 
of the diagram. Secondly, we distributed 28 points along 
the isotherm at e/ksT — 600 for pcr"^ within [1 — 8], and 
tested the mechanical stability of all Bravais crystals and 
the three non-Bravais structures considered by running a 
short NVT simulation for each specific structure. In par- 
ticular, for lattices with free parameters, typically simu- 
lations at 5 different values for each free parameter were 
carried out. This large set of simulations reduced sig- 
nificantly the number of crystal structures for which we 
subsequently computed free energies accurately: in par- 
ticular, at pa^ = 4 we calculated free energies for the SC, 
BCC, FCC, R and H. At pa^ = 5, for the SC, H, and R. 
At pcr3 ^ 7^ foj. ^Yie ST, BCT, BCO, H, and R. 

To determine the thermodynamically-stable crystals, 
we obtained free energies for the set of mechanically- 
stable candidate struct ures. We employed standard ther- 
modynamic integratioiP^'^Sl using the Einstein crystal as 
a reference state and a 20-point Gauss-Legendre numer- 
ical integration. To assess accuracy, particularly at high 
densities, we carried out four integrations for each crys- 
tal structure, at values of aa^ /ksT equal to 50, 100, 150 
and 200, where a is the Einstein-crystal spring constant. 
The calculated free energies are thus averages over four 
integrations. Our NVT simulations typically consisted of 
500-700 particles, the actual number depending on the 
crystal structure. In particular, we used 500 spheres for 
the FCC, 686 for the BCC, 512 for the H and SC, and 
648 for the BCT. Occasionally, simulations with about 
1000 particles were run at a few points distributed in the 
T — p space, but we did not find any significant finite-size 
effects. 

The determination of coexistence points involved NPT 
simulations to compute the local equation of state P{p) 
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3 3 

pa Pa I z 

FIG. 2: Temperature-density and temperature-pressure cuts of the phase diagram of Hertzian spheres. Errors are smaller 
than the size of the circles. Solid lines are a guide to the eye. We performed accurate free-energy calculations for pa^ < 7. 
Broken lines thus indicate approximate phase boundaries. The inset on the left zooms in the region around the F-FCC-BCC 
triple point, and the inset on the right compares the Fourier transform of the Hertz (solid line) and 3D overlap (dotted line) 
potentials. See section |V] for details. 



TABLE V: Some phase-coexistence points. Standard deviations for densities are below 5 10 



Phases 1-2 


e 


10^ ksT/e 








Phases 1-2 


e 


10^ fcsT/e 




Pl(T^ 




F - FCC 


12000 


8.33310^2 


12.34 


1.000 


1 


104 


F - H 


280 


3.571 


497.4 


3.644 


3.709 


F - FCC 


1000 


1.000 


14.74 


1.145 


1 


252 


F - H 


250 


4.000 


474.5 


3.764 


3.811 


F - FCC 


400 


2.500 


17.78 


1.295 


1 


389 


F - H 


228 


4.386 


501.2 


4.046 


4.050 


F - FCC 


200 


5.000 


24.15 


1.531 


1 


602 


H - F 


270 


3.704 


767.4 


4.579 


4.518 


F - BCC 


125 


8.000 


44.75 


2.003 


2 


028 


H - F 


240 


4.167 


617.1 


4.361 


4.325 


F - BCC 


113.1 


8.842 


71.17 


2.398 


2 


398 


H - SC 


12000 


8.33310-2 


38290 


4.725 


4.822 


FCC - BCC 


12000 


8.33310"^ 


5912 


2.214 


2 


279 


H - SC 


1000 


1.000 


3132 


4.693 


4.785 


FCC - BCC 


1000 


1.000 


494.6 


2.216 


2 


275 


H - SC 


400 


2.500 


1189 


4.597 


4.676 


FCC - BCC 


400 


2.500 


194.3 


2.206 


2 


252 


H - SC 


270 


3.704 


768.6 


4.518 


4.586 


FCC - BCC 


200 


5.000 


88.35 


2.148 


2 


175 


SC - F 


270 


3.704 


802.4 


4.674 


4.677 


FCC - BCC 


150 


6.667 


58.23 


2.071 


2 


086 


SC - F 


325 


3.077 


1215 


5.184 


5.227 


FCC - BCC 


130 


7.692 


41.28 


1.959 


1 


966 


SC - F 


400 


2.500 


1647 


5.419 


5.478 


BCC - F 


280 


3.571 


491.9 


3.540 


3 


625 


SC - BCT 


12000 


8.33310-2 


48510 


5.379 


5.506 


BCC - F 


200 


5.000 


310.1 


3.355 


3 


427 


SC - BCT 


1000 


1.000 


4091 


5.402 


5.526 


BCC - F 


140 


7.143 


164.6 


2.997 


3 


045 


F - BCT 


400 


2.500 


1665 


5.508 


5.558 


BCC - F 


125 


8.000 


123.9 


2.808 


2 


840 


F - BCT 


350 


2.857 


1595 


5.760 


5.799 


BCC - H 


12000 


8.33310-2 


20470 


3.453 


3 


686 


F - BCT 


310 


3.226 


1574 


6.076 


6.105 


BCC - H 


1000 


1.000 


1724 


3.477 


3 


701 


BCT - F 


325 


3.077 


2069 


6.750 


6.796 


BCC - H 


400 


2.500 


699.5 


3.516 


3 


699 


BCT - F 


300 


3.333 


1683 


6.381 


6.384 



around each pair of guesses for the equilibrium densities. 
We used 30 simulation runs for each guess and fitted a 10- 
degree polynomial function to the data. The free energy 
of the solid phase at the estimated density was calculated 
via integration from the Einstein crystal as explained in 
the paragraph above. We computed the free energy of 
the fluid by integrating energies from 800-particle NVT 



simulations along a constant density path starting at the 
isotherm at ksT/e = 1/105. This isotherm consisted 
of polynomial fits to the equation of state of the fluid 
for eight different equally-spaced intervals of pa'^ within 
the range [0 — 8]. We obtained free energies along the 
isotherm by integrating the equation of state from the 
ideal-gas limit at zero density. We generated as many 
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simulation points as necessary to keep standard devia- 
tions of integrals below 0.01%. Accuracy was absolutely 
important also when comparing free energies of compet- 
ing crystals, since differences were often below 5 lO^^fcsT 
per particle. 

Once the local equations of state P{p) and chemical 
potentials fi{p) at fixed T were determined, we searched 
for the densities pi and p2 that satisfy Pi{pi) = P2{p2) 
and pi{pi) — /i2(p2)- We then checked that pi and p2 fall 
within the interval for which P{p) and p{p) were fitted. 
Otherwise, we used pi and p2 as new guesses for the coex- 
istence densities, and repeated the procedure explained 
in the previous paragraph. 



IV. PHASE DIAGRAM 

According to the criterion established by Likos in 
ref. [m soft potentials that are purely repulsive and 
bounded give rise to re-entrant melting when its Fourier 
transform exists and it does not have negative compo- 
nents. If the Fourier transform oscillates around zero, 
cluster solids appear instead. The Fourier transform of 
the Hertz potential is definite positive and, indeed, our 
results confirm the re-entrant melting scenario, as Fig. [2] 
shows. Table |V] summarizes phase coexistence data. 

The phase diagram of Hertzian spheres shows that the 
fluid phase freezes upon compression to form an FCC 
crystal which, at higher densities, turns into a BCC struc- 
ture. However, the BCC packing is favored over the FCC 
at high temperatures because particles in the former have 
a higher vibrational entropy. An inset in Fig. [2] shows 
in detail the area around the F-FCC-BCC triple point, 
which occurs at ksT/e = 7.6910-3. The FCC phase 
does not have a maximum freezing point for the Hertz 
potential, as opposed to the GCM.^ The fluid re-enters 
at densities larger than the maximum freezing point at 
fcsT/e = 8.84 lO'^ and pa^ = 2.40. At lower temper- 
atures and with increasing density, we find multiple re- 
entrant melting and four more crystals: hexagonal, sim- 
ple cubic, body-centered-tetragonal and trigonal. Simple 
cubic crystals of one component are rare, because they 
are usually mechanically unstable. 

With regard to the H and BCT structures, the geom- 
etry of their primitive cells is not unique, as explained 
in section |II A[ Both crystals have one degree of free- 
dom, namely the aspect ratio of the primitive cell c/a 
(see Table One expects that the optimal value of c/a, 
at which the free energy reaches a minimum, is in gen- 
eral density dependent. We show this in Fig. [3] where 
we plotted free energies per particle as a function of the 
aspect ratio. For the H crystal, the minimum in free en- 
ergy occurs at c/a « 0.84, independently of temperature. 
Within the range of densities at which the H structure 
is stable, we do not observe any significant density de- 
pendence either, as the free-energy curve is relatively fiat 
around the minimum. However, this is not the case for 
the BCT. As shown in the inset of Fig. [3] the maximum 




0.82 0.84 0.86 0.760 0.765 
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FIG. 3: Free energies per particle of the hexagonal (left) 
and body-centered tetragonal (right) crystals, as a function 
of the respective aspect ratio c/a of the primitive unit cell. 
Free-energy sets have been displaced arbitrarily by /o for each 
specific point in the phase diagram, as indicated in the leg- 
ends. Errors are smaller than the size of the symbols. The 
optimal c/a, at which the free energy is at a minimum, is 
density dependent for the BCT structure. The inset shows 
that the density behavior of the optimal c/a has inflection 
points such that aspect ratios close to that of the BCT at the 
maximum melting point are favored. Lines are a guide to the 
eye. 

freezing point, occurring at pa^ — 6.35, leaves its "im- 
print" on the density dependence of the optimal aspect 
ratio of the BCT primitive cell, as there is a clustering of 
aspect ratios close to c/a — 0.762, which is the value at 
the maximum freezing point. 



V. DISCUSSION 

Very recently, Prestipino et al^ calculated the 
ground-state structures for the Hertz potential. They 
found that the non-bravais cI16 and A5 lattices minimize 
the energy at densities within the interval [4.3 — 6.4]. We 
did not include the cII6 and A5 in the set of candidate 
structures, and hence we cannot tell whether the cI16 
and A5 remain stable at T > 0. However, we stress 
that entropic contributions are not negligible even at the 
lowest temperatures that we studied. For example, at 
pa^ = 5, the trigonal crystal with a/w = 0.46 yields 
the lowest energy of all the candidate structures consid- 
ered in this work. Yet, the simple cubic structure has 
the lowest free energy; it is thus stabilized by entropy. 
Furthermore, as discussed in detail below, the phase dia- 
gram and both structural and dynamical anomalies seen 
in the fiuid are interconnected. For example, the maxi- 
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mum melting point of the BCT lies at pa^ = 6.35, which 
is approximately the density at which the radial distri- 
bution function and the diffusion coefficient reach local 
extrema. But at T = a transition between the A5 and 
the BCT happens to be close to that same densitji^H (the 
A5 becomes unstable at pa^ — 6.386), which makes an 
extremum occurring at the same location in the phase 
diagram unlikely. The fact that there is such a strong 
correlation between the shape of the computed phase di- 
agram and the dynamics of the fluid provides indirect 
support for the assumption that we have indeed identi- 
fied the stable crystal phases at T > 0. 

In the next paragraphs we discuss a few interesting as- 
pects of the fluid phase. The re-entrance of the fluid at 
high densities and low temperatures (e.g. at pa^ — 7.5 
and kBT/e = 2.5 10"'^; see Fig. [2| appears to be different 
in nature from that of a hard-sphere fluid, which has a 
rugged potential-energy landscape and yet is a fluid be- 
cause there is an extended and connected region in config- 
uration space where the potential energy is strictly zero 
(and therefore flat). For the Hertz model at the afore- 
mentioned density, the typical inter-particle distance in 
the liquid phase is of the order of 0.5a and the potential is 
at less than 20% of its maximum value; also, at the afore- 
mentioned temperature, the relative radial displacement 
of two particles by as little as 2.5 10""^ cr will change the 
potential energy of such a pair by an amount that is com- 
parable to kgT. Hence, the potential-energy landscape is 
far from being flat at these conditions. We suggest that 
the reason why Hertzian spheres (and other potentials 
showing re-entrant melting, like the GCM) melt at high 
densities and low temperatures has to do with the fact 
that the Fourier transform of the potential goes to zero 
rapidly. Indeed, SiitcP^ (see also ref.dSJ has shown rigor- 
ously that integrable bounded potentials whose Fourier 
transform is non-negative and vanishes above a wave 
number Kq have an infinite number of continuously de- 
generate ground states above a well-defined threshold 
density. This means that, for this class of pair poten- 
tials, the potential energy is completely flat along certain 
directions in configuration space. Although the Hertz 
potential does not satisfy the Siito criteria (its Fourier 
transform does not vanish above any Kq), its Fourier 
transform does decay cx k~'^ for wave vectors larger than 
^ lOfj^^ (see inset in Fig. |2|. We speculate that, as a 
result, certain collective motions in the dense phase cost 
very little potential energy. It is the softness of the pair 
potential that causes melting into a rather peculiar liquid 
for which particle motions should be strongly correlated. 

In Fig. |4] we show other aspects of the peculiarity of 
the dynamics of the fluid phase. The diffusion coefficient 
D, which decreases monotonically with increasing den- 
sity in simple liquids, shows a clear non-monotonic be- 
havior in the Hertzian fluid. Moreover, the minimum in 
D at kBT/e — 10^^ coincides with the maximum freez- 
ing point at pa^ — 2 AO, and a maximum in the first 
peak of the g{r) (see inset in Fig. |4]). Also the maxima 
in D for the three temperatures shown in Fig. [4] are lo- 
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FIG. 4: Diffusion coefficient of the fluid phase obtained by 
molecular dynamics as a function of density at three temper- 
ature values. Breaks in the lines joining the simulation data 
indicate that points in between were identified as crystalline 
states. The inset shows the radial distribution function g{r) 
for the isotherm at ksT/e — 0.01 and the following values of 
pa^: 0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, 8, 10, 12. The highest 
peak in each g{r) is marked with a black dot. In the set of 
black dots, the local maximum and minimum occur approxi- 
mately at the densities at which the diffusion coefficient has 
a minimum (p(T^=2.40) and a maximum {pa'^— 6.35), respec- 
tively. 



cated at the same density as the maximum freezing point 
at pa'^ — 6.35 and a minimum in the first peak of g{r). 
These multiple structural and dynamical anomalies, and 
their interrelation, has been seen before in bounded po- 
tentials,'22122l in liquids with water-like anomalie^^ and 
colloids with short-ranged attractions.!^ This suggests 
that crystallization and fluid dynamics are inextricably 
connected, and that this connection is not independent of 
the underlying crystal structure. For instance, the pres- 
ence of the H and SC crystals does not imprint peaks in 
the diffusion coefficient as the BCT does. We think that 
this is related to the fact that single-particle mobility 
is affected by differences in packing efficiency: the in- 
trinsic volume of the BCT primitive cell is much smaller 
than that of the BCC, H and SC structures. Besides, 
even though particle motions in the low-temperature liq- 
uid phase are strongly correlated, the liquid does not ap- 
pear to be particularly glassy. In fact, within the range 
of densities studied, we observe spontaneous crystalliza- 
tion rather than structural arrest upon cooling even at 
temperatures below ksT < 10""*. Further study of the 
dynamics of the high-density, supercooled liquid would 
therefore be most interesting. 
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Another interesting feature of the system at high den- 
sities and low temperatures is its freezing behavior. We 
have used free-energy calculations to trace freezing curves 
up to a density of pa^ ~ 7. Even beyond that density, 
there appears to be a succession of other freezing transi- 
tions. Although we have not traced the accurate phase 
boundaries completely for pa^ > 7, we have calculated 
free energies at pa^ = 9 and ksT/e = 10^'^, and found 
the stable structure to be the BCT with c/a ~ 1. Ad- 
ditionally, we have found that crystallization occurs at 
low temperatures and at densities up to pa^ — 12. In 
other words: as far as we can tell, the low-temperature 
phase is always crystalline, but the structure of the sta- 
ble crystalline phase changes as the density is increased. 
There is also evidence^ at T = supporting this obser- 
vation. In fact, the ground-state behavior is interesting 
in the context of a result reported by Torquato and Still- 
inger.-^ on the basis of duality relations between a soft 
potential and its Fourier transform, these authors found 
a one-dimensional system, the overlap potential, that ex- 
hibits an infinite number of phase transitions between 
periodic ground states over the entire density range. One 
can then wonder whether the Hertz potential also shows 
an unbounded number of phase transitions between peri- 
odic structures. Interestingly, the Hertz potential and the 
overlap potential in ID and 3D, Vmir) = e(l — r/a) and 
Vsoir) = e(l-(3/2)(r/a) + (l/2)(r/cr)3) for r < a, stud- 
ied in ref. 33, belong to the same class of non-negative. 



bounded functions, and both have oscillatory, decaying 
Fourier transforms (see ref. l34l and inset in Fig.[2|. 

In summary, we have shown that the Hertz potential 
gives rise to a phase diagram with multiple re-entrant 
melting transitions and a succession of Bravais crystals 
spanning a wide range of densities. In addition, the dif- 
fusion coefficient of the re-entrant fluid shows unusual 
non-monotonic curves with increasing density. This rich 
behavior together with the simplicity of the potential's 
functional form makes the Hertz model an attractive one 
for the study of kinetic phenomena in soft-core systems, 
like martcnsitic transformations and supercooled dynam- 
ics. 
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